/* Include files */
#include "math.h"
#include "mex.h"

/* Include to sort : http://stackoverflow.com/questions/1787996/c-library-function-to-do-sort */
/*
#include <stdio.h>
#include <stdlib.h>
*/

/* Checks whether inputs are matrices */
#define IS_REAL_2D_FULL_DOUBLE(P) (!mxIsComplex(P) && \
mxGetNumberOfDimensions(P) == 2 && !mxIsSparse(P) && mxIsDouble(P))
#define IS_REAL_SCALAR(P) (IS_REAL_2D_FULL_DOUBLE(P) && mxGetNumberOfElements(P) == 1)

/*************************************************/
/* MAIN                                          */
/*************************************************/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* Define output matrix */
	#define Share    	plhs[0]
	#define ShareSt 	plhs[1]

	/* Define input matrices */
	#define Firm 		prhs[0]
	#define Market      prhs[1]
	#define J 			prhs[2]

	int i, j, m;

	/* Check the number of inputs is correct (Aoffer,Baccept,QbarA) */
	if(nrhs < 3)
		mexErrMsgTxt("Too few input arguments.");
	else if(nrhs > 3)
		mexErrMsgTxt("Too many input arguments.");

	if(!IS_REAL_2D_FULL_DOUBLE(Firm))
		mexErrMsgTxt("Firm must be real 2D full double array.");
	if(!IS_REAL_2D_FULL_DOUBLE(Market))
		mexErrMsgTxt("Market must be real 2D full double array.");
	if(!IS_REAL_2D_FULL_DOUBLE(J))
		mexErrMsgTxt("J must be real 2D full double array.");

	/* Dimensions and pointers for inputs */
	int N = mxGetM(Firm);
	int M = mxGetM(Market);

	double *Firm_pt 	= mxGetPr(Firm);
	double *Market_pt 	= mxGetPr(Market);
	double *J_pt 		= mxGetPr(J);

	int Jm = J_pt[0]*M;
	int Jn = J_pt[0]*N;

	/* Output matrix */
	Share = mxCreateDoubleMatrix(Jm,1,mxREAL);
	double *Share_pt = mxGetPr(Share);
	for (j = 0; j < Jm; j++)
	{
		Share_pt[j] = 0;
	}

	ShareSt = mxCreateDoubleMatrix(Jn,1,mxREAL);
	double *ShareSt_pt = mxGetPr(ShareSt);
	for (i = 0; i < N; i++)
	{
		ShareSt_pt[i] = 0;
	}

/***************************************************************/
/* Begin Filling       */
/***************************************************************/
/*
Begin, End : market index begins and ends for parcels
Foreach match in FirmMatch : add 1 if j,m
Assumes that each firm can enter each market (J_pt[0] count of firms)
*/

	int shareID = 0;
	int shareStID = 0;
	int beg = 0;
	int end = 0;
	int Jind = 0;
	int denom = 0;

	for(m = 0; m < M; m++)
	{
		beg = Market_pt[m + M*1];
		end = Market_pt[m + M*2];
		for (j = beg - 1; j < end; j++)
		{			
			if (Firm_pt[j] > 0)
			{
				shareID = Firm_pt[j] + m*J_pt[0] - 1;
				denom = Market_pt[m + M*0];
				Share_pt[shareID] = Share_pt[shareID] + 1/(float)denom;
			}
		}
	}
	for(m = 0; m < M; m++)
	{
		i = 0; Jind = 0;
		beg = Market_pt[m + M*3];
		end = Market_pt[m + M*4];
		for (j = beg - 1; j < end; j++)
		{			
			shareID = i + m*J_pt[0];  		/* Firms are stacked (1:J)' */
			shareStID = (beg - 1) + Jind*J_pt[0] + i;
			ShareSt_pt[shareStID] = Share_pt[shareID];
			i++;
			if (i == J_pt[0])
			{
				i = 0;
				Jind++;
			}
		}		
	}
}
